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ABSTRACT 

Mass outflows driven by stars and active galactic nuclei are a key element in 
many current models of galaxy evolution. They may produce the observed black 
hole-galaxy mass relation and regulate and quench both star formation in the 
host galaxy and black hole accretion. However, observational evidence of such 
feedback processes through outflows of the bulk of the star forming molecular gas 
is still scarce. Here we report the detection of massive molecular outflows, traced 
by the hydroxyl molecule (OH), in far- infrared spectra of ULIRGs obtained with 
Herschel-PACS as part of the SHINING key project. In some of these objects 
the (terminal) outflow velocities exceed 1000 km/s, and their outflow rates (up 
to ~1200 M0/yr) are several times larger than their star formation rates. We 
compare the outflow signatures in different types of ULIRGs and in starburst 
galaxies to address the issue of the energy source (AGN or starburst) of these 
outflows. We report preliminary evidence that ULIRGs with a higher AGN 
luminosity (and higher AGN contribution to Ljfi) have higher terminal velocities 
and shorter gas depletion time scales. The outflows in the observed ULIRGs 
are able to expel the cold gas reservoirs from the centres of these objects within 
^10*^-10^ years. 



^Max-Planck- Institute for Extraterrestrial Physics (MPE), GiessenbachstraBe 1, 85748 Garching, Ger- 
many 

^Universidad de Alcala de Henares, 28871 Alcala de Henares, Madrid, Spain 

■^Department of Astronomy, University of Maryland, College Park, MD 20742, USA 

^Naval Research Laboratory, Remote Sensing Division, 4555 Overlook Ave SW, Washington, DC 20375, 
USA 

^Tel Aviv University, Sackler School of Physics & Astronomy, Ramat Aviv 69978, Israel 
f^Oxford University, Dept. of Astrophysics, Oxford 0X1 3RH, UK 

^INAF-Osservatorio astronomico di Roma, via Frascati 33, 00040 Monteporzio Catone, Italy 



- 2 - 



Subject headings: ISM: jets and outflows — galaxies: active — galaxies: evolution 
— galaxies: starburst 



1. Introduction 

Gas-rich galaxy merging may trigger major starbursts, lead to the formation of elliptical 
galaxies, and account for the formation and growth of supermassive black holes (BHs; e.g., 
Sanders et al. 1988; Hopkins et al. 2009). This merger-driven evolutionary scenario starts 
with a completely obscured ultraluminous infrared galaxy (ULIRG). As the system evolves, 
the obscuring gas and dust is gradually dispersed, giving rise to dusty QSOs and finally to 
completely exposed QSOs. Powerful winds, driven by the central quasar or the surrounding 
starburst, have been invoked to quench the growth of both the BH and spheroidal component 
and explain the tight BH-spheroid mass relation (e.g., di Matteo et al. 2005, Murray et al. 
2005). These winds are purported to inhibit star formation in the merger remnants ( "negative 
mechanical feedback"; e.g. Veilleux et al. 2005 for a review), and to create a population of 
red gas-poor ellipticals, thereby explaining the bimodal color distribution observed in large 
galaxy surveys (e.g., Kauffmann et al. 2003). 

Finding observational evidence of such feedback processes in action is one of the main 
challenges of current cxtragalactic astronomy. While outflows have been observed frequently 
in many starbursts and QSOs, so far they have been detected mostly in the ionized and 
neutral atomic gas component. To inhibit star formation in the host galaxy, outflows have 
to affect the molecular gas out of which stars form. Few detections of molecular outflows 
have been reported so far (e.g. Baan et al. 1989, Walter et al. 2002, Sakamoto et al. 
2009). In this Letter we demonstrate that far- infrared molecular spectroscopy with Herschel- 
PACS of (ultra-) luminous infrared galaxies is providing a breakthrough in identifying and 
analyzing massive molecular outflows. Our recent OH-absorption observations have revealed 
a >1000 km/s molecular outflow in the closest quasar known, Mrk231 (Fischer et al. 2010). 
Independent, spatially resolved CO-emission observations of Mrk231 with the IRAM/PdB 
mm- wave interferometer (Feruglio et al. 2010) have conflrmed this outflow with inferred 
mass outflow rates of ~300-2200 M0/yr, signiflcantly larger than the current star formation 
rate (SFR^lOOMg/yr) in the host galaxy. We now show that molecular outflows are indeed 
a common phenomenon in many of the luminous major mergers in our sample, reaching 
outflow velocities of 1000 km/s and outflow rates up to ~1000 M0/yr in some of them. 



^Herschel is an ESA space observatory with science instruments provided by European- led Principal 
Investigator consortia and with important participation from NASA. 
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2. Observations and Data Reduction 

The data presented here are part of the Herschel guaranteed time key program SHIN- 
ING, a study of the far- infrared properties of the ISM in starbursts, Seyfert galaxies and 
infrared luminous galaxies. Here we present velocity-resolved line profiles of the OH 79/im 
cross ladder, ground-state doublet for a small subset of our ULIRG sample and a comparison 
starburst galaxy. The observations were obtained with the PACS far-infrared spectrometer 
(Poglitsch et al. 2010) on board Herschel (Pilbratt et al. 2010). Our data set also contains 
other OH transitions, but the 79/im line was observed first and is common to all objects, 
as it can be obtained simultaneously (in second order) with the [CII] 157/im line (in first 
order). At 79/im the PACS resolution is ~140 km/s. The data reduction was done using the 
standard PACS reduction and calibration pipeline (ipipe) included in HIPE 5.0. However, 
for the final calibration we normalized the spectra to the telescope flux (which dominates the 
total signal, except for NGC253) and re-calibrated it with a reference telescope spectrum 
obtained from dedicated Neptune observations during the Herschel performance verification 
phase. All of our objects (except NGC253) are point sources for PACS. In the following we 
use the spectrum of the central 9"x9" spatial pixel (spaxel) only, applying the point source 
correction factors (PSF losses) as given in the PACS documentation. We have verified this 
approach by comparing the resulting continuum fiux density level to the continuum level of 
all 25 spaxels combined (which is free of PSF losses and pointing uncertainties). In all cases 
the agreement is excellent, however the central spaxel alone provides better S/N. We note for 
completeness that for NGC253 the total OH 79/im line profile summed over all 5x5 spaxels 
yields emission, consistent with the ISO-LWS observations by Bradford et al. (1998). 

In a next step we have performed a continuum (spline) fit. Due to the limited wavelength 
coverage these fits are somewhat subjective. To help define continuum points and potential 
additional spectral features (such as the H2O absorption line at 78.74/im, indicated with an 
arrow in Fig. [1]) we have used our full range spectra of Arp220 and NGC4418. These two 
sources will be analyzed in detail in forthcoming papers, but preliminary data points for 
Arp220 are included in Figures 2 and 3. We note here that NGC4418 shows signatures of 
an infiow. 



3. Targets 

For this first study of outfiow signatures in our data we use a sub-sample that is mainly 
constrained by the observing schedule of Herschel, but that covers a broad range of AGN and 
starburst activity, including a starburst template (NGC253), a cold, starburst dominated 
ULIRG (IRAS 17208-0014,) warm ULIRGs (S25/S60>0.1) and/or ULIRGs with strong AGN 
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Fig. 1. — Observed PACS spectra (continuum normalized) of the OH transition at 79 /xm 
(grey). Overplotted are the low velocity (dotted) and high velocity (dashed) fit components 
and the total fit (full). The arrow indicates the rest position of H2O 423-812 . The dash-dotted 
line for IRAS 14378 shows the observed spectrum of the OH transition at 119 fim for this 
object. 



Table 1: Target properties, outflow rates and outflow velocities (la uncertainties in brackets) 



Source 


SFR 

Me/yr 


a"" 

% 


lOi^Lo 


gas 

lO^Mo 


Me/yr 


y%ak 

km/s 


^85% 

km/s 


vf 

max 

km/s 


Mrk231 


101 (15) 


71 (11) 


28 (4) 


4.2 (1.3) 


1 1 qn+4700 


-600 


-660 


-1170 


IRAS 08572+3915 


42 (6) 


72 (11) 


12 (2) 


1.3 (0.4) 


970l?^r 


-700 


-740 


-1260 


IRAS 13120-5453 


168 (25) 


9 (1.4) 


1.8 (0.3) 


5.8 (1.7) 


130+^f 


-520 


-600 


-860 


IRAS 14378-3651 


>79 


<45 


<7.2 


4.2 (1.3) 


7/in+2200 
'^'-'-550 


-800 


-860 


-1170 


IRAS 17208-0014 


274 (41) 


11 (1.7) 


3.4 (0.5) 


12.2 (3.7) 


qn+270 


-100 


-170 


-370 


NGC253 


1.7 (0.3) 








0.7 (0.2) 


1 a+i.S 


-75 


-130 


-280 



^ Fraction of the AGN contribution to Lboi, where Li,oi = l.lSxL//?; ^ Gas mass (taken from Gracia- 
Carpio et al. 2011); ^ mass outflow rate (see footnote of Table [2]); peak velocity of the blue shifted 
high velocity component (relative to systemic velocities); velocity for which 85% of the outflowing 
gas has lower (absolute) velocities; ^ Terminal velocity. Estimated uncertainty for all velocities: 
±150km/s. 
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contributions (Mrk231, IRAS 13120-5453, IRAS 14378-3651), and a heavily obscured ULIRG 
(IRAS 08572-1-3915), which hosts a powerful AGN (e.g. Veilleux et al. 2009, V09). 

4. Results and Discussion 

Figure [T] shows the (continuum normalized) OH 79/im line spectra for all objects. For 
NGC 253 we show the central spatial pixel only. The Mrk231 spectrum is taken from Fischer 
et al. (2010) and is repeated here for completeness. In all cases we detect P-Cygni profiles 
typical of outflows, with blue shifted absorption and red shifted emission features, of more 
than 1000 km/s in some cases. For comparison with the literature we list in Table [1] various 
measures of outflow velocities (relative to the system velocity of the blue component of the 
OH doublet): the peak velocity (vpeak), the maximum (terminal) velocity (ymax), and V85% 
for which 85% of the outflowing gas has lower (absolute) velocities. The uncertainties of 
these velocities are dominated by the uncertainties in the continuum definition, the S/N 
ratio in the spectra, and the spectral resolution. We estimate an overall error of ±150 km/s. 
The S/N in the IRAS 14378-3651 spectrum is relatively low, but the high terminal velocity 
is confirmed by the OH 119 fim transition (see Figured]). 

These high OH outflow velocities may be the long-sought conclusive evidence of powerful 
mechanical feedback from vigorous star formation and/or accreting central black holes. The 
possible feedback and outflow mechanisms (e.g. winds from SNe, radiation pressure) are 
debated in the literature (see, e.g., the review by Veilleux et al. 2005). It is not clear, 
if such mechanisms could indeed be sufficient to power outflows that are strong enough 
to significantly affect the host galaxy and to actually quench the star formation in these 
objects. It is also unclear from the models whether it is possible to distinguish AGN driven 
outflows from stellar-driven outflows observationally (see, e.g., Hopkins & Elvis 2010). In 
the following we adopt an empirical approach with our new data. 

4.1. Are the strong outflows we observe driven by the AGN rather than by 

the star formation in these objects? 

Rupke et al. (2005a,b,c) and Krug, Rupke & Veilleux (2010) have studied large samples 
of AGN and star forming galaxies in neutral gas (blueshifted optical Na ID 5890, 5896 A 
absorption features). They found that, for fixed SFR, ULIRGs with higher AGN fractions 
have higher neutral gas outflow velocitites, reaching velocities well above 1000 km/s in some 
broad-line AGN (see also, e.g., Heckman et al. 2000, Martin 2005, 2006, Thacker et al. 
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2006). Theoretical models predict that SNe driven outflows cannot reach velocities higher 
than 500-600 km/s (e.g. Martin 2005, Thacker et all. 2006). Predictions of v^ax from 
models of outflows driven by radiation pressure from a starburst (or AGN) are, however, 
less certain. The terminal velocity we measure in the OH outflow of our starburst template 
NGC253 is ~300 km/s. IRAS 17208-0014, a starburst dominated ULIRG with little AGN 
contribution has only a slightly higher Vmax (370 km/s). In significant contrast to this the 
terminal velocities of the OH outflow in the two AGN-dominated ULIRGs (Mrk231 and 
IRAS 08572+3915) are well above 1000 km/s. Thus the OH outflow velocity could be a 
very promising tool to distinguish AGN driven outflows from starburst driven outflows, with 
AGN-dominated outflows reaching much higher velocities. 

In Figure[2]we compare the terminal outflow velocities to the star formation rates (SFRs) 
and AGN-luminosities of our objects. SFRs are calculated from the IR luminosities (based 
on the calibration of Kennicutt 1998, but using a Chabrier IMF rather than a Salpeter 
IMF), i.e. SFR = (1-a) x lO^^^L/ij, applying AGN correction factors a (the fraction of the 
contribution from the AGN to L^o/, where l^hoi = l-15xL/ij as in V09). These factors are 
individually derived from Spitzer mid-IR diagnostics provided by the QUEST programme 
(V09). For IRAS 13120-5453 and IRAS 14378-3651 (which are not part of V09) we calculated 
AGN fractions using method 1 from V09, based on the observed [OIV]/[NeII] ratio (and 
the observed upper limit for IRAS 14378-3651) from Farrah et al. (2007). AGN luminosities 
are then calculated as Lagn = ctxLhoi- AGN luminosities and SFRs are listed in Tabled] 

For star formation driven outflows one might expect the outflow velocities to scale with 
the star formation rates (e.g. Tremonti et al. 2007). We do not see such a correlation 
(Figure [21 upper panel). Instead, we see a rough correlation of Vmax with Lagn (Figure El 
lower panel), consistent with the idea of the high velocity outflows being powered mostly by 
(radiation pressure from) the AGN. 

Very energetic outflows have been found in radio galaxies at high z (e.g. McCarthy et 
al. 1996, Best et al. 2000, Nesvadba et al. 2009). Radio jets could in principle be driving 
the outflows we see. However, nearly all of the objects presented here are radio quiet, so 
this energy source can safely be assumed to be negligible. Only in Mrk231 does a radio jet 
contribute to the NalD outflow (Rupke & Veilleux 2011). 
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Table 2: (Preliminary) model fit results (errors are discussed in the text) 



Source 


component 


TDC 


n(OH)f„ 


r 




<N(OH)>g 






pc 


lO'^^cm"^ 




degrees 


lO^^cm-^ 


Mrk231 


HVC^ 


105 


7 


0.6 


90 


3 




LVC^ 


115 


5 


1.0 


61 


5 


IRAS 08572+3915 


HVC 


110 


20 


0.2 


77 


4 




LVC 


110 


40 


0.3 


69 


4 


IRAS 13120-5453 


HVC 


210 


9 


0.3 


43 


3 




LVC 


210 


1 


0.5 


90 


2 


IRAS 14378-3651 


HVC 


100 


3 


0.5 


90 


1 




LVC 


120 


2 


1.0 


90 


1 


IRAS 17208-0014 


LVC 


110 


6 


1.0 


57 


2 


NGC253 


LVC 


90 


0.5 


1.0 


51 


1 



^ high velocity component; ^ low velocity component; inner shell radius; OH density at Rm; 
filling factor; ^ half opening angle (i.e. 9 = 90° for a full Att coverage); ^ OH column density N(OH)x/. 
These values are used to compute the mass outflow rate in Table [1] 

M - Mgas/tdyn 47r X n{OH)^n / x{OH) X X Rl^x f X g XV, where tdyn ^ R/v, g is a function 
of the opening angle 9 and x{OH) = 5 x 10^^ is the OH abundance relative to H2. Outflow rates in 
Tabled] are the sum of the components (HVC+LVC). A more detailed description of a refined model 
will be given in a future paper. 
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4.2. Does the outflow carry suflicient molecular gas to remove the star 
formation fuel and actualy quench the star formation? 

To compute mass outflow rates (M) for comparison with SFRs we have modeled the ob- 
served spectra using the radiative transfer code described in Gonzalez-Alfonso & Cernicharo 
(1999). The outflow is modeled as concentric expanding shells around a nuclear continuum 
source, allowing for each source one to three components with different velocity gradients and 
distances to the central source. Such geometry reproduces quite naturally the red-shifted 
emission (which is produced in the receding cocoons), and the blue-shifted absorption (from 
the approaching parts). 

The density proflle for each velocity component is determined through mass conservation 
assuming a stationary outflow. With the exception of NGC253, where we flt the nuclear 
continuum emission from the observed PACS spectrum, the nuclear FIR in all other sources 
is modeled as the warm component fltted in Mrk231 (Gonzalez- Alfonso et al. 2010). Free 
parameters are the inner and outer radii and velocity fleld of each velocity component, 
the OH density at the inner radius, the covering factor of the continuum FIR source, and 
the solid angle of the outflow. Because the OH energy levels are radiatively pumped in 
the outflows, transitions at different wavelengths and energy levels, in combination with 
continuum component flts, yield crucial information about the radial location where the lines 
are formed. Besides the 79 /xm doublet, the 119 (also ground-state) and 65 fim (-Eiower = 290 
K) doublets were also observed in all our sources, except for IRAS 17208. The excitation 
of the high-lying 65 fim doublet requires high far-IR radiation densities, indicating that the 
line is tracing the inner, highly excited, molecular region of the outflowing material, close 
to the nuclear source (Gonzalez- Alfonso et al. 2008). Models that account for the 65 /im 
line absorption simultaneously reproduce the two ground-state lines as well, indicating that 
the three observed OH lines - in contrast to the atomic NalD lines, which are extended on 
galactic scales in some cases - are mostly sensitive to the inner ~ 0.5 kpc of the outflow. 
For IRAS 17208 only the 79 /im line is currently available, but we have assumed a similarly 
compact outflow in OH. The OH columns are typically several 10^^ cm~^. In Table [2] we 
list the flt results. Figure [1] shows the various components of the spectral flts on top of the 
observed spectra. 

We derive mass-loss rates M (see Table [T] and footnote of Table [2]) from the OH density 
profile of the outflowing gas, the gas velocity in each component, and the adopted OH 
abundance. The resulting mass- loss rates are higher than 700 Mo/yr in those sources with 
high Vmax- The uncertainty in M is dominated by the adopted OH abundance, the properties 
of the underlying FIR continuum source, and the outflow geometry. We adopt an overall 
uncertainty of 4, derived from a study of the dependency of the flt results on the various 
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model assumptions. The OH abundance relative to H2 is assumed to be 5 x 10~^, based on 
modeling of multi-transition OH observations of the Galactic giant molecular cloud SgrB2 
(Goicoechea & Cernicharo 2002). While OH abundances could reach values significantly 
lower in different environments (e.g. in some of the models of Sternberg & Dalgarno 1995 
and Meijerink & Spaans 2005, or the observations by Watson et al. 1985), we expect a 
relatively high OH abundance close to the nuclear region (Gonzalez- Alfonso et al. 2008) and 
note that M increases with decreasing OH abundance. On the other hand, even a high OH 
abundance of 10~^ leads to mass loss rates of several hundred M-q/ji. Further details on the 
modeling will be given in a forthcoming paper. 

The starburst-dominated objects of our sample (NGC253, IRAS 13120-5453, and IRAS 
17208-0014) have values of M similar to their star formation rates. In contrast to this, 
the AGN-dominated and/or warm ULIRGs, as well as the heavily obscured ULIRG, have 
outflow rates that are factors ~4-20 larger than their SFR (see Tabled] and Figure [3], upper 
panel). The outflows we observe in the ULIRGs - if continued at the current rate - are able 
to completely expel the cold gas reservoirs from the centres of these objects within ~10^- 
10^ years (Figure [31 lower panel). The gas depletion time is roughly inversely correlated 
with the AGN luminosity (and AGN fraction). This is consistent with the merger scenario 
where the highest outflow rates are a short-lived, late, AGN-dominated stage in the merger 
evolution. We note, however, that the lowest time scales we find (few x 10^ years) appear to 
be inconsistent with star formation/galaxy evolution models. We may have underestimated 
the time scales: a higher gas mass by a factor 2 and a higher OH abundance by a factor 
2 would result in 4 times longer time scales (which would still be within the error bars of 
figure 3). A fundamental change in the basic model assumptions, like the geometry of the 
outflow, could also change the outflow rates. Our current models are consistent with wide 
opening angles for the outflowing gas, but this is hard to prove with PACS spectra alone (see 
below). For Mrk231, however, the high outflow rate is consistent (within a factor 2) with 
the independent CO observations (PdBI mm-interferometry) by Feruglio et al. 2010. Very 
recently the high rates and the wide angle of the outflow have also been confirmed by optical 
spectroscopy through the analysis of broad wings of Ha and NalD absorption (Rupke & 
Veilleux 2011). Finally, the outflow may not continue at such a high rate, since - as gas 
and dust are removed - the outflows may be quenched as the radiation pressure decreases. 
Therefore, our depletion time scales are time scales rather then absolute durations. 

We do not have the spatial resolution with Herschel-PACS to study the spatial distribu- 
tion of the outflowing material directly (although observations of additional OH transitions 
with higher energy level, e.g. at 84 and 163 /im, could help to further constrain the geometry 
of the models). True spatial information on the molecular component can only come from 
interferometric mm-observations (e.g. Feruglio et al. 2010). While the observed spectra 
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Fig. 2. — Maximum outflow velocities (terminal velocities) as function of star formation 
rate (upper panel), and AGN luminosity (lower panel). The asterisk denotes NGC253, the 
triangle Arp 220. 
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Fig. 3. — upper panel: the ratio of the mass outflow rate to the star formation rate versus 
SFR; lower panel: depletion time scale versus AGN luminosity. Symbols as in Figure El 
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shown here clearly exhibit high velocities, powerful molecular outflows, independent of our 
modelling, the conclusions based on outflow rates and AGN fractions need further confir- 
mation from increased statistics. Our ongoing observations will cover a much larger sample 
of both starburst and AGN dominated (U)LlRGs than shown here. Gombined with mm- 
interferometric follow up observations we will be able to better constrain our OH outflow 
models and further investigate the surprisingly low depletion time scales for some of the 
objects. The final sample should then allow us to study potential trends in the outflow 
characteristics. Examples of such trends could be a different or tighter correlation of the 
outflow velocity with AGN luminosity than with star formation rate, or the extent to which 
AGN fractions and outflow velocities reflect different merger stages with evolving outflows 
creating lower and lower covering factors to the AGN. 
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